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It is likely that ultra-high energy cosmic rays contain a significant component of heavy or in- 
termediate mass nuclei. The propagation of ultra-high energy nuclei through cosmic radiation 
backgrounds is more complicated than that of protons and its study has required the use of Monte 
Carlo techniques. We present an analytic method for calculating the spectrum and the composition 
qq \ at Earth of ultra-high energy cosmic rays which start out as heavy nuclei from their extragalactic 

, sources. The results obtained are in good agreement with those obtained using numerical methods. 
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Recent observations by the High Resolution Fly's Eye PJ and the Pierre Auger Observatory (PAO) [2] have estab- 
lished that at energies above ~6x 10 19 eV, the spectrum of ultra-high energy cosmic rays (UHECR) is suppressed 
significantly below the power law form that holds at lower energies. At such high energies, protons are expected to 
,—— / interact catastrophically with photons in the cosmic microwave background (CMB) with an attenuation length of 
^ ! 0(100) Mpc - the Greisen-Zatsepin-Kuzmin (GZK) 'horizon' [IHH]. 

However the observed attenuation does not 
Qh' necessarily imply that the UHECR are protons since heavy nuclei too will lose energy through photodisintegration by 
q , photons of the cosmic infrared background (CIB) with a similar energy loss length [E0]- 

The possibility that UHECRs contain a substantial fraction of heavy or intermediate mass nuclei is interesting to 
consider for two reasons. Firstly, recent measurements of UHECR shower profiles at the PAO |8j reveal a gradual 
decrease in the the average depth of shower maximum above ~ 2 x 10 18 eV, indicating increasing dominance by heavy 
nuclei at the highest energies. Secondly, the maximum energy to which cosmic rays can be accelerated in ther sources 
is proportional to their electric charge, making it less challenging for plausible astrophysical accelerators to produce 
£> ' heavy nuclei at the highest observed energies Q . 

Recently, the PAO collaboration has also found a correlation between the arrival directions of UHECRs with energies 
C*"} ■ greater than 6 x 10 19 eV and nearby active galactic nuclei within ~ 75 Mpc This does suggest that most of 

these particles are not deflected by galactic or extragalactic magnetic fields by more than a few degrees. As cosmic 
ray nuclei have higher electric charge than protons and consequently suffer bigger deflections by magnetic fields, it 
is tempting to use this observation to argue in favor of a proton-dominated UHECR spectrum. UHECR nuclei can 
however be significantly disintegrated during propagation, leading to a population of lighter particles observed at 
Earth, thus reducing the impact of magnetic fields (especially in the Galaxy) on UHECR deflection. Furthermore, 
there are considerable uncertainties in the magnitude and structure of galactic and extragalactic magnetic fields. 
Hence it is still pertinent to consider the possibility that the primary cosmic rays are heavy or intermediate mass 
nuclei and study their propagation through radiation backgrounds, in order to quantify the expected composition and 
spectrum at Earth. 

After a l ong hiatus, there has been quite a bit of work in the past decade on the propagation of UHECR nuclei [ll|, 

m m ei ei ee d m, ei he m, si m, hi m m, but we propose here a different a pp r ° acn to thc problem. 

While the complex process of the photodisintegration of UHECR nuclei into lighter nuclei and nucleons has so far 
been addressed using Monte Carlo techniques, we develop an analytic description of this phenomenon. Our primary 
motivation is to reduce the computation time needed to calculate the observed UHECR spectrum and composition 
for a given injected spectrum, composition and cosmic ray source distribution. 

The rest of this article is organized as follows. In Sec. [H] we outline the coupled differential equations governing the 
population of nuclear species during the photodisintegration process and find an analytic solution through physically 
justified simplifying assumptions. In Sec. IIIII we apply this to a cascade initiated by ultra-high energy iron nuclei, 
obtaining expressions for the spatial distribution of each of the species populated in the cascade. In Sec. HVi we 
compare the results of our analytic approach to those obtained using Monte Carlo methods and demonstrate their 
close agreement. In Sec. [V] we present our conclusions. 
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II. THE ANALYTIC FORMULATION 



As UHECR nuclei propagate through intergalactic space, they interact with cosmic radiation backgrounds, frag- 
menting into lighter nuclei and nucleons at a rate: 

A* m y /"» den(e) r 2Ee / Am P c2 
RA,z,t p , tn = 2E2 g2 J aee<T A ,z,i p ,in( £ )> (!) 

where A and Z are the atomic number and charge of the nucleus, i p and i n are the numbers of protons and neutrons 
broken off from a nucleus in the interaction, n(e) is the density of background photons of energy e, and <JA,z,i p ,i n (e') 
is the appropriate cross section (see Appendix • 

An exact analytic treatment of the propagation of UHECR nuclei would need to take into account all of the many 
possible decay chains (ie. all values of i p and i n for each nuclear species), as we have done previously using a Monte 
Carlo technique [HI, EH, H^j. The resulting differential equations are non-trivial to solve, and there would be no 
significant computational advantage over the Monte Carlo approach. Denoting the number of nuclei with atomic 
number A by Na, the differential equation describing the population of this state is given by: 

dN A N A N A _ N A+1 N A+2 

+ i + ; + . . . — + +..., (Z) 

ClL Ia^A-1 'A-»A-2 lA+l^A lA+2^A 

where L is the distance traveled and is the interaction length for the state i to disintegrate into state j (For a 
description of how the interaction lengths, Zj_ >j, are calculated, see Appendix [A"f . To begin with, we consider the 
simplified case in which only single nucleon loss processes occur. This reduces the number of states in the system 
dramatically, along with the number of possible transitions. Eq. ([2|) then simplifies to: 

dN A , N A N A+1 ^ 



dL l A I 



A+l 

where a slight change of notation has been introduced — what was written earlier as Ia->a-i is now denoted simply 
as I a (since a particle in state A+l may decay only into state A in our simple model). 

The solution to this set of coupled differential equations, constrained by the initial conditions N n (L = 0) ^ and 
Na{L = 0) = (for 4/n), is given by (see Appendix IB"]) : 

^=£mt-'» p (-^) If r+r- m 

m=A x ' p=A{^m) 

where the interaction lengths, I a, denote those for single nucleon loss {A — ► A — 1). 

The degree of agreement of this result with that obtained by Monte Carlo can be improved greatly by redefining 
I a in terms of an effective interaction length: 

which accounts for multi-nucleon loss processes as well. 

From Eq. Q , the average number of various nuclear species can be obtained as a function of the distance travelled 
by the parent cosmic ray. At the high energies of interest, we can sensibly ignore energy losses due to pair creation 
(negligible relative to photodisintegration) and the cosmological redshift (negligible for the nearby sources which 
dominate) . Hence we can simply relate the energy of the nuclei to the energy E n of their parent particle (of mass 
number n): Ea = E n (A/n). From now on we will denote the initial particle energy simply as E so the fraction of 
particles in state A with energy Ea{= EA/n) after the parent particle (of mass number n) has propagated a distance 
L, is given by 



An analogous equation has been used to describe near-threshold pion production by ultra-high energy protons prop- 
agating through the cosmic microwave background (CMB) [281 ] . 
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FIG. 1: The average surviving fraction of various nuclei (left panel) and the average number of nucleons disassociated from 
each primary nucleus (right panel) as a function of the distance propagated from an extragalactic source of 10 20 eV iron nuclei. 



III. THE NUMBER AND SPATIAL DISTRIBUTION OF SECONDARY NUCLEI AND PROTONS 

In each photo-disintegration interaction, a proton is produced (neutrons are produced too but decay quickly into 
protons). The differential equation governing the population of these protons, denoted by Ni, is: 

dNijE) = N n (L,E n ) JV n _i(L,£; n -i) N 2 (L,E 2 ) 

dL L n (E n ) L„_ 1 (£^_i) L 2 {E 2 ) ' [ ' 

where E — E n /n. Note that all protons produced during a cascade initiated by a primary cosmic ray of mass number 
n with energy E have an energy E / A n . Their energy distribution is however subsequently altered through interactions 
with the cosmic background photons. 

Integrating Eq. (fTj) over L, we obtain an expression for N±, at distance V and energy E\(= E/A n ): 



N 1 (L',E 1 )= / diV 

Jo „._<; 



\N m (L,E m ) 



m=2 v ' 



(8) 



From this expression it is clear that over sufficiently large distances, each of the different species contributes equally 
to the proton injection spectrum, providing n protons in total. Following from Eq. the distribution of each species 
at a given distance and energy is given by: 

\ ' m=q x v ' ' p=q(^m) v ' f\ v/ 

where E m = E(A m /A n ) and E p = E(A p /A n ). 

As an illustration we apply this to the specific case of 10 20 eV iron nuclei, adopting the same CIB spectrum as in 
our previous work [16l |. The surviving fraction of various nuclear species as a function of distance from the source is 
shown in the left frame of Fig. [TJ while the right frame shows the average number of nucleons broken off the injected 
nuclei. From such distribution functions, the average composition of cosmic rays arriving at Earth from any set of 
sources with a specified spatial distribution and energy spectrum can easily be derived. 

Focussing now on the distribution of all particles with a given energy, E, the distribution of such particles is just 
the sum of the terms in Eq. §§§ over all species (requiring different energy parent nuclei at source): 

dNML,E) _" N q (L,E q ) 

dL ^N n (0,E n ) ' 1 ' 

q—l 

where now E q = E(A n /A q ) (note the inversion). 

If the particles are injected over a range of distances L, from to £ max , with a distribution described by the 
normalized function d(L), then the number of a particular species, q, at an energy, E, is given by: 

N q (E q )= [ Lm ™ dL M^ d[L) , (11) 



N n (0,E n ) 
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FIG. 2: The fraction of particles which start out as 10 20 eV iron nuclei in representative states (A=50, 40, 30, and 20) as a 
function of the distance propagated, calculated using both analytic and Monte Carlo techniques. 

If particles are injected with a differential energy spectrum dN/dE cx E~ a , the total number of all particles with 
energy E q is just 



^(E q )^±^N q (E q )=±(^ 



dE 

q=l 

where I is the lightest species considered in the cascade. 



q=l 



N q (E), 



(12) 



IV. COMPARISON OF THE ANALYTIC SOLUTION WITH MONTE CARLO RESULTS 



With a knowledge of the effective interaction lengths of each state, La(E), the population of any species of interest 
after propagating over a specified distance can be determined following Eq. ([5]). We will now compare the results of 
this approach with those from the Monte Carlo calculation described in our previous work [l6| • 

In Fig. [21 we show the population of five representative nuclear species, as a function of distance from the source, 
as calculated using both analytic and Monte Carlo techniques. We find excellent agreement, especially for heavier 
nuclear species. The small degree of disagreement for lighter nuclei seen in Fig. [2] is the result of multi-nucleon loss 
processes, the effect of which is most significant over long decay chains. 

We continue the comparison of the analytic and Monte Carlo methods in Fig. [3l where we plot the number of 
protons produced through photodisintegration as a function of the distance propagated. The results from the two 
methods are virtually indistinguishable. 

We can also use our analytic method to calculate the average composition of the UHECR spectrum after propagation. 
The average mass number, (A), of the particles arriving with energy, E, is given by, 



(ME q )) 



Y,UA q N q {E q ) 
N tot (E q ) ■ 



(13) 



Integrating the spatial distribution function over a given set of cosmic ray sources, the total number of each species 
in the cosmic ray spectrum can be obtained. A calculation of the number of proton secondaries from photodisintegra- 
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FIG. 3: The number of protons produced through the photodisintegration of 10 eV iron nuclei, as a function of the distance 
propagated, calculated using both analytic and Monte Carlo techniques. 
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FIG. 4: The ultra-high energy cosmic ray spectrum (left) and average composition (right) calculated using both analytic and 
Monte Carlo techniques, for the case of iron nuclei injected by homogeneously distributed sources with dN/dE oc E~ 2 up to a 
maximum energy of 10 22 eV. The black curves show the MC results with pair creation losses and redshift effects included. 

tion processes which arrive at Earth is slightly more involved, and requires consideration of both the production of 
secondary protons through photodisintegration and their subsequent pair creation and pion production energy losses 
on cosmic photon backgrounds. We can estimate the required quantity for protons: 

fL P (E 1 ) 

Ni{E\) = / dL'Ni(L', Ei), (14) 

Jo 

where L p is the energy loss length for protons. 

In Fig. [H we compare the spectrum and composition after propagation found using both our analytic technique 
and the Monte Carlo programme, for the case of iron nuclei injected by a homogeneous distribution of sources with 
a spectrum dN/dE oc E~ 2 up to a maximum energy of 10 22 eV. Due to the energy range of interest in these plots, 
only sources up to a redshift z=l were considered. The agreement between the two techniques is rather good — the 
analytic method captures all of the essential features in the results from the Monte Carlo (in which pair creation losses 
and redshift effects were also included). 1 It is seen that the simplifying approximations made in our analytic approach, 
namely neglecting pair creation and redshift energy losses, are quite appropriate for energies above ~ 10 19 eV. 



1 The decrease in (A) between 10 20 3 and 10 21,2 eV can be understood as follows. Since the Fe nuclei arc injected to a maximum energy 
of 10 22 eV, the protons produced by photo-disintegration have a maximum energy 56 times smaller, i. e. 10 20 - 25 eV. Above this energy, 
only heavier particles may contribute. As shown in Fig. 6 of our previous paper |23|| . at 10 20 ' 3 eV the interaction lengths of low mass 
nuclei have decreased to a minimum already while the interaction lengths of high mass nuclei have just started decreasing. Hence as 
the energy is increased even further, the contribution of the heavier nuclei to the arriving flux decreases, hence so does (A). 
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V. CONCLUSION 



We have presented an analytic method to compute the UHECR energy spectrum and chemical composition at 
Earth resulting from the propagation through cosmic radiation backgrounds of heavy nuclei injected by extragalactic 
sources. This reduces dramatically the computation time compared to the usually used Monte Carlo programme 
and provides insight into the physical reasons for the results obtained. These results can be confronted with the 
observational data, to enable the identity of the primary particles to be established. 
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APPENDIX A: INTERACTION LENGTHS 



To carry out the calculations described, first the interaction lengths for the photodisintegration of nuclei need to 
be calculated. These are given by: 

jl_ _ r de^e) (Al) 



l A ^A-i 2£ 2 

where A is the atomic number of the nucleus, i is the number of nucleons broken off from a nucleus in the interaction, 
n(e) is the differential number density of background photons of energy e, and a Ai i(e') is the appropriate cross section. 
We model the photodisintegration cross sections with the parameterization of Refs. [H, 0, EH : 

^E d W i - 1 e- 2 ( £ - £ -) 2 / A 'e + (eti„-)e_(e 1 ), e thr < e < e 1; i = 1, 2 
CA,i(e) = { C s d6+(ei)B_(e max )/(e max - ei), a < e < e max (A2) 

0; £ ^ ^max 

where £j, C, e p,i an d A^ are parameters whose values are obtained by fitting to nuclear data and the integrated 
cross-section is 



E d ee / a(e) de = 2vr 2 - % {A = 60 {A "/^ mb-MeV, (A3) 

' 4-7reo m p c z A A 



with Z being the charge of the nucleus. The function Wi is given by: 

A - ei 



TP, 



erf S l + erf 



Ai/V2 J V A 4 /V2 



(A4) 



Here, 6 + (x) and 8_(x) are the Heaviside step functions, t\ — 30 MeV, e max = 150 MeV^and the threshold energy 
for a given process is in most cases is ethr ~ i x 10 MeV (values are tabulated in Ref. [13|). These cross sections 
are dominated by the giant dipole resonance which peaks in the energy range ~ 10 — 30 MeV; at higher energies, 
quasi-deuteron emission is the main process. 

Nuclei with A = 5 — 9 are very unstable and quickly decay to helium and lighter elements. We approximate this 
behaviour by setting Ia^4 = for A = 5 — 9. 



APPENDIX B: DEVELOPMENT OF A SOLUTION 

Consider the differential equation, 



dN A N A N A+1 

CH L A l A+ i 
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This may be written as, 



exp 



I \ d 



I A dl 



exp [ j^ n a 



N A 



+1 



Ia 



+1 



so that 



Ni = exp [ -— ) I dl exp ( — 



I \ N A+1 



with the initial conditions of the system being, N n (L = 0) ^ and Na{L = 0) = (for A^n) 
Let us assume that the solution is: 



N A (l) 
N n (0) 



m—A 



I, 



i - 1 ' 



where n is the initial state the particles arc in. The following will be a 'proof by induction'. 
If the above is true for A, then it is true also for A + 1 : 



N A+ i(l) 
N n (0) 



— E ^A+llm 



A-1 



exp 



m=A+l 



I 



£ n 



i 



p=A+l(^m) 



1 - / ' 



which using Eq. (|B3|) gives 



A^(Q 

iV„(0) 



= «p(-l)/di«p(l) e c 



Ia 



i-A-2 



exp 



E i 

m=A+l 
n 

E w 



A-2 



n-A-1 



m=A+l 

-1 



J__ _L, rxp(- — 
t^4 t r , 



n 



n 

1 



J p=A+l&m) 



, exp ( - — 



exp 



m=A+l 



r n 



i 



cexp i - — 



Applying the boundary condition Na(L = 0) = 0, we get 

n n 

c = e wr 1 - 1 n 

m=A+l p=A(=£m) 

For Eq. (|B4[1 and Eq. (|B6j) to be equivalent, it is necessary that, 

n n 1 

e^- 1 n 7z 

m—A 



p-A(^m) 



(B2) 



(B3) 



(B4) 



(B5) 



(B6) 



(B7) 



(B8) 



This final expression may be recognised as an expression of the 'Vandermond determinant' with the final column 
repeated. Since this is indeed true, our assumed solution |B4I must also be correct. 
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